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Recently, lancu and Triantafyllopoulos have proposed a hierarchy of evolution equations in QCD 
at high energy which generalises previous approaches by including the effects on gluon number 
fluctuations and thus the pomeron loops. In this paper, we present the first numerical simulations 
of the Langevin equation which reproduces that hierarchy. This equation is formally the Balitsky- 
Kovchegov equation supplemented with a noise term accounting for the relevant fluctuations. In 
agreement with theoretical predictions, we find that the effects of the fluctuations is to reduce the 
saturation exponent and to induce geometric scaling violations at high energy. 
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. I. INTRODUCTION 

^' 

^ • The question of the high-energy limit of QCD is a longstanding problem. The resummation of the leading logarithmic 
' ' contributions at high-energy, well-known as the BFKL equation Jj, predicts a power growth s'^ of the scattering 
t . amplitude T. In the large- iVc approximation, the resummation of these leading logarithmic contributions is equivalent 
to the evolution of a wavefunction considered as a system of qq dipoles. 
, At sufficiently high energy, this violates the Froissart bound and the unitarisation effects have to be taken into 
• account. In the large-A'c approximation, it has been shown that the relevant way to introduce these unitarisation 
[ effects is to include multiple-scattering effects in the collision between external dipoles and evolved target. This leads 
I to an infinite hierarchy of equations 0, taking into account linear BFKL growth and saturation effects in the 

■ high-density regime. In the mean-field approximation, valid e.g. for dense targets such as heavy nuclei, the infinite 
. hierarchy resumes to a single non- linear equation known as the Balitsky-Kovchegov (BK) equation ^|^. 

Jlf^ ' Recently, one has noticed 0,11] that, for dilute systems, one also have to include the effects of fluctuations of the 
number of gluons. It has been subsequently recognised, by lancu and Triantafyllopoulos that these contributions 
f| are not included in the original Balitsky/JIMWLK hierarchy but can be added to it. The additional terms describe 
O-i gluon splittings in the target wavefunction and act as a seed for the pomeron loops. The resulting infinite set of 
ft) • equations is equivalent, up to a coarse-graining approximation, to a stochastic Langevin equation which formally have 
Qj ' the form of the BK equation supplemented with a noise term responsible for fluctuations 9]. 

<~i In the saddle-point approximation, Munier and Peschanski have shown ^3 that the BK equation is equivalent to the 

L| Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) equation [H well-known in statistical physics. As a consequence 
• i-H of this equivalence, one knows that the asymptotic solutions of the BK equation takes the form of a traveling wave 

■ and leads to geometric scaling 0,0] i.e. to the fact that the scattering amplitude depends only on k"^ IQl{Y) where 
^ , is the transverse momentum of the dipole and Qs{Y) is the rapidity-dependent saturation scale. If we includes the 

■ - - ' noise term, the equation becomes equivalent to the stochastic FKPP (sFKPP) equation. In this type of equation, 
geometric scaling is observed event by event but the averaged amplitude leads to violations of geometric scaling as 
the events are diffused around their average position. 

The equivalence between QCD and statistical-physics equations is actually an approximation. If we go beyond 
the saddle point approximation, i.e. if we look at the exact QCD evolution equation, we get additional non-local 
contributions. These can be studied using the same techniques, leading to slightly modifled asymptotic properties 
and to a different infrared behaviour. In addition, the analytic results H, Q obtained in statistical physics are only 
valid for asymptotically small values of the strong coupling constant as . Only a numerical analysis can give relevant 
results in the physical domain. So far the only numerical studies of the QCD equations including fluctuations |l3l | has 
been restricted to models without transverse coordinate dependence. 

The aim of this paper is to study numerically the properties of the solutions of the Langevin equation describing 
the QCD amplitude and to compare with the results obtained from BK evolution. We shall check, in the case of the 
QCD Langevin equation and for real-life values of the physical parameters, that the effects of the fluctuations are well 
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(a) linear term (b) fan diagram (c) fluctuation 

FIG. 1: These figures show the contribution to the evolution of T'^': a linear growth proportional to T'^', a merging term 
proportional to T*-^^ and a splitting term proportional to T. 



predicted by the sFKPP-inspired approximation 9] at large rapidities and in the tail of the front. Basically, we are 
going to study the effect of fluctuations on the dispersion of the events, the Y dependence of the saturation scale and 
the form of the wavefront. Let us emphasise that our studies are directly related to the QCD equation, hence we do 
not solve the sFKPP equation or any other model such as the effect of a cut-off [3| on the amplitude. 

The structure of the paper is as follow. We shall first explain the basic considerations leading to the hierarchy 
of evolution equations including fluctuations and introduce the Langevin equation. We shall then then expose the 
results on the asymptotic solutions obtained from the equivalence with statistical physics equations. We then explain 
how we solve numerically the QCD equation and the results obtained from this analysis. We finally conclude and give 
the perspectives of this paper. 



II. FLUCTUATIONS IN QCD EVOLUTION 



In this section, we briefly recall the basis steps yielding to the Langevin equation for the rapidity evolution of QCD 
amplitudes. The limit we want to consider is large Y = log(s), and as small enough to allow for a perturbative 
treatment. We shall work in the framework of the large- iVc approximation and introduce d = ^"^^ . In order to study 
the rapidity evolution of a target, we shall probe it with a set of projectile made of qq dipoles. 

It is now well-known 0, [l5j | that, due to nonlinear unitarisation effects, it is not sufficient to consider interaction 
with one single dipole under the form of a BFKL ladder. In other words, it is not sufficient to consider a simple two- 
gluon exchange in the t channel and the associated amplitude T, but one has to take into account the higher-order 
correlations (T(xi, Y) . . . T(x„, Y)). These nonlinear effects introduces 0, |^ in the evolution of the n dipole system 
a contribution from the n -\- 1 dipole sector. This translates the fact that unitarity effects in dense systems take the 
form of recombinations. 

It has recently been shown 0i S Pi that in addition to the unitarity corrections relevant in the dense regime, it is 
also of prime importance to take into account gluon number fluctuation terms if one wants an acceptable description 
of the dilute systems. If one put all these effects together, lancu and Triantafyllopoulos obtain an infinite hierarchy 
of equations for all dipole correlations. If T^*^^ (xi , yi ; . . . ; x^ , ; V) is the correlation for the interaction with k dipoles 
of coordinates (xi,yi), (xfe,yfe), this hierarchy takes the form 

9yTW(xi,yi;...;xfc,yfc;y) 
= d( d?zY. ( ^""'w^'^' ^2 [r^'^(xi,yi;...;x„z;...;Xfc,yfc;y)+rW(xi,yi;...;z,y,;...;Xfc,yfc;y) 

-T(^) (xi, yi; . . . ; Xfe, yfe; Y) - T^^^+i) (xi, yi; . . . ; x„ z; z, y, ; . . . ; x^, y^; 
In X: ' y 7/^'+^^' T(^-^)(x„ yi; . . . ; x„ y,+,; . . . ; x., y.; F) ^(^)(y, - x,+0, (1) 



where k runs over all positive integers. These equations are obtained if one makes the approximation that the 
amplitude is proportional to ag(x— y)^n(x, y) where n is the gluon number density. The proportionality constants 
appears as the unknown factor k in the fluctuation term. This approximation is correct in the BFKL or in the dilute 
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regime. The constant k is expected to be of order 1. In fact, this approximation has been relaxed in subsequent 
pubhcations (see second paper in Ref. [3 and HI) but the non-locahty of the ensuing equation is too comphcated 
to be accounted for in numerical simulations. This is why throughout this paper we shall restrict ourselves to the 
simplified hierarchy (Q. In fact, in what follows we shall consider a simplified version of equation obtained after 
a coarse-graining in impact-parameter space. 

To make things clearer, let us discuss the evolution equation for the correlations of two dipoles T^^^ 



9yT(2) (xi,yi;x2,y2;r) 

(xi -yi' 



(xi - z)2(z - yi)2 



, - 2 (Xl 

-I- aa„K 



-yi)^(x2 -y2)^ 

(xi - y2)2 



T(2)(xi,yi;x2,z;y)+r(2)(xi,yi;z,y2;r)-T(2)(xi,yi;X2,y2;r) 
-T(3)(xi,yi;x2,z;z,y2;r) + (1^2) 
r(i)(xi,y2;r)5(2)(yi-x2). 



(2) 



We distinguish three different types of contributions to the rapidity evolution of T^^^^. The first one is the BFKL linear 
term, proportional to T^^^^ itself. It corresponds to the diagram of figure ^a) and contributes to the elaboration of 
a BFKL ladder. Then, the term proportional to T'^^\ shown in figure ^b), resums fan diagrams, responsible for the 
saturation effects. Finally, the last term (figure ^c)) corresponds to the fluctuations and introduces a contribution 
from T in the evolution of T^'^\ 

It is interesting to check the regimes in which these terms are expected to be important. This is easily done by 
considering the power of as associated with each graph. Since T represents at least a two-gluon exchange, each power 
of T carries 4 vertices i.e. a factor a^. One thus easily see that the terms proportional to T^^^^ or T'-^'> both carry 
an additional factor as while the fluctuation term, proportional to T has an extra factor a^. This means that the 
merging contribution becomes important when T^'') ^ T ^ 1. It takes into account saturation effects in the dense 
regime. Moreover, the ffuctuation term becomes important when a^T ~ T^^-*, i.e. when T ~ a^. One thus expects 
large effects due to the fluctuations in the dilute regime. As an extreme case, if one starts with a single dipole, we 
have T^'^^ = for fc > 1 and the fluctuation term is needed to obtain a nonzero value for the correlations. 

Note that in the hierarchy, the physical roles of the QCD parameters d and are very different. While d 
always multiplies Y and hence sets the scale for the rapidity evolution, a'^ measures the strength of the elementary 
dipole-dipole interaction and thus sets the scale for the fluctuations. In fact, in the numerical simulations, we found 
convenient to flx the value of the QCD coupling constant (fixing simultaneously the values of d and a^) and to use 
the freedom in the choice of the parameter k as a method to vary the strength of the dipole interaction. 

If we neglect the fluctuation term in this equation, one recovers the large- TVc version of the JIMWLK hierarchy 
0, which is coherent with the fact that it is then expected to describe dense systems such a heavy nuclei. It is also 
interesting to notice that, while the JIMWLK equation at large resums just fan diagrams, the complete hierarchy 
including the fluctuations also allows for pomeron loops. As a simple example, in two steps of evolution, T^-^"^ can 
fluctuate into T'^^ and merge back into T^^^. 

If we perform a coarse-graining approximation in impact parameter, the fluctuation term becomes local and we can 
perform the Fourier transform to momentum space (see Ref. for technical details). By applying the same methods 
as those exposed in appendix the resulting hierarchy of equation in momentum space can be rewritten under the 
form of a Langevin equation 



dYT{k,Y) = d 



dp2 

/o 

where the noise i^{k, Y) satisfies 



p'^T{p,Y) - k^T{k,Y) k'^T{k,Y) 



|p2-/c2| 



dT^{k) + dyj2KalT{k)v{k, F), 



(3) 



(i/(fc,y))==0 and {v{k,Y)v{k' ,Y')) = — 5{Y - Y') k 5{k - k'). 

air 

This equation has to be understood as the continuum limit with the Ito prescription which means that if we perform 
rapidity steps of size e, labelled by an index j , it can be rewritten 



Tl+l!^ — ZiM ^aK® Tj{k) - dTf{k) + d^2KalTj{k)vj{k), 



(4) 



where K ®T denotes the convolution with the linear kernel and the noise Vi satisfies 



{vj{k)) = and {vj{k)vf{k')) = —SjfkS{k - k'). 
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III. TRAVELING WAVES AND ASYMPTOTIC PROPERTIES 



Now that we have presented the equation describing high-energy evolution of QCD amphtude, let us quote what 
are the expected properties of its solutions. 

The linear term in (|SJ) can be rewritten x(— with L = log(fe/fco) and x{l) = 2?/'(l) — ^'(7) — V'll — 7) is 
the eigenvalues of the BFKL kernel. If one expand this operator to second order around 7 = 1/2, the equation © 
becomes after a linear change of variables and a rescaling of T 



dtu{x, t) = d'^u{x, t) + u{x, t) — u^{x, t) + yj2rju{x, t)v{x, t). 



(5) 



This equation is the stochastic Fisher-Kolmogorov-Petrovsky-Piscounov (sFKPP) equation [ij well known in statis- 
tical physics. This analogy helps to understand most of the asymptotic properties of the evolution equation |M 0| . 

Let us first consider the equation without noise (FKPP equation) , equivalent to the BK equation in QCD [l^l • This 
equation is known to admits travelling waves as asymptotic solutions [13 . In the framework of the BK equation, this 
means that at high rapidity, the solution takes the form of a travelling front at high k 



T{k,Y) 



Qs{Y) 



(6) 



where 



QliY) 



(7) 



is the saturation scale. Both the critical slope 7c ~ 0.6275 and the critical velocity Vc ~ 4.8836q; can be determined 
from the BFKL kernel x{l)- Iii addition, the form of the traveling front looses traces of the initial condition^. 

If we now switch on the noise term this picture gets modified and one has to consider the sFKPP equation. First 
of all, if one takes one single event of the Langevin equation, the asymptotic average velocity of the propagating front 
is smaller than Vr uM 



V ^ Vc 



21og2(a2K) ' 



(8) 



As indicated in the above equation, this analytic result is only valid for asymptotically small values of the noise 
strength a^K. The shape of the front is also modified by the noise. It has even been shown |0 that, for the sFKPP 
equation, if we start with an initial condition which is always equal to 1 for x < Xq and to for x > xi^ then this 
compacity property is preserved by evolution. The extension of the tail is then proportional to log(K). More precisely, 
the tail of the front behaves like 



A\ log(a^K)| sin 



Ctt 



log 



|log(a2^)| ^\Ql{Y) 



exp 



'Ic log 



If we now consider a set of events, the effect of the noise will have a diffusive effect. This means that at rapidity Y , 
the events will show similar fronts distributed normally around Qs{Y) = {Qs{Y)). The width of this dispersion is 



Aiog[g^(r)]« v^disay 



(9) 



typical of Gaussian white noises. It has been numerically found 14] that, for asymptotically small a^K, the front 
diffusion coefficient scales like 



1 



diff ' 



\\og\a^,K)\ 



The point at this stage is that, although event by event geometric scaling is satisfied, if we average over multiple 
events, it is destroyed as we typically get an average amplitude behaving like a front traveling at a constant speed 
but spreading like \/Y Q . 



^ This is only true for steep enough initial conditions which, for the case of QCD, is ensured by colour transparency. 
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In the case of the equation (O, if one assume a front compact and saturating to a constant value, lancu and 
Triantafyllopoulos have shown that the average value of the amplitude is 




where z = log(fc^/(5^) is the position relative to the wavefront, ~ DdisaY is the diffusion of the events in- 
duced by fluctuations and erfc(a;) is the complementary error function. The normalisation factor N is given fixed 
to |x"(l/2)\/l + 8x(l/2)/x"(l/2) ~ 6.98 according to the transformation leading to the sFKPP equation. Finally, 
their model also predicts that, at asymptotic rapidities, the tail of the front satisfies 

(T) ~ (r^) ~ • • • ~ (T") ~ . . . (11) 

It is important to note that our purpose is not only to test these predictions. The predictions (|HJl and Q, obtained 
from the equivalence with the sFKPP equation, are only valid in the limit where the noise is very small i.e. in the limit 
of small and unphysical (or, equivalently, small k). In this work, we shall study the solutions of the QCD equation 
for physical values of the strong coupling. In addition, by analysing the complete QCD equation, without using the 
saddle point approximation, we shall also take into account the fact that, in the infrared region, the amplitude exhibits 
a logarithmic behaviour instead of saturating to a constant. 



IV. NUMERICAL SOLUTIONS 
A. Description of the method 

The aim of this paper is to solve numerically the evolution equation Q without simplifying it to the sFKPP 
equation and without using a model such as a inserting cut-off. To achieve this task we need two ingredients: an 
efficient way to deal with the linear-kernel convolution, and an efficient way to take into account the noise term. 

In order to observe the front propagation to large values of k as energy increases, we shall work with the variable 
L — \og{k'^ /kg) where ko is a fixed soft reference scale. The integration in ||3J) is then performed by discretising space 
between Lmin and L^ax- The compactness of the front in the presence of a noise term ensures that we can neglect 
the ultraviolet contribution coming from T{L > Lmax)- This is a crucial property which accelerates the convergence 
of the numerical algorithm w.r.t. the case of the BK equation. In the infrared region, one knows that the amplitude 
shows a logarithmic divergence i.e. T{k) « log(l/fc) -I- est for small k. Assuming that behaviour for L < Lmin, one 
can compute analytically the contribution to the linear term coming from the infrared. 

The treatment of the noise term is a little bit more tricky. The main problem is that, due to the presence of 
the square root in the pre-factor, one has to ensure that the noise does not produces negative values of T. The 
method we have used is the one proposed by Levine and Pechenik | 18i |. The technique is to consider the probability 
associated with the stochastic process. For the one under consideration, the probability that, in a rapidity interval 
Y, the amplitude goes from Tq to T satisfies the following equation 

dyPiY, T- To) = vd^ [TP{Y, T; To)] , 

with rj — Kaa^/ {nV AL) (Ai is the step for discretising the space dimension) and the initial condition 

F(o,r;To) = ,5(r-ro). 

One can solve this differential equation and then impose that the amplitude remains positive. The final result is [T^ 

Once we know this probability, one can define the cumulative probability 

Fto;Ay{T)= [ dT' P{AY,r;To), (13) 
No- 
where the prescription on the inferior integration bound is such that the first term in (|12|) is included. 
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The method to apply the noise in then straightforward: if u represent a uniform random variable on the interval 
[0, 1], T = -Fjjj^Ay(^) distributed according to the probability (|12|l . We can thus define an intermediate amplitude 

-^noisc(^) "^T(fc);Ay t^^^-^-l ' 

where u{k) is the uniformly-distributed random variable at position k. Once we have taken into account the fluctuation 
term we apply the remaining part of the evolution equation in order to obtain the evolved amplitude 

T(r + Ay, k) = T„oisc(fc) + Ay [aK ® T^oiscik) ~ aT^^isM] ■ 

Finally, the integration H13|) cannot be performed analytically. In order to simplify the numerical treatment of the 
noise term which already involves heavy computations, we have approximated the T > part of the probability H12() 
by a suitable Gaussian function leading to an error function for the cumulative frequency H13|) . 

In our numerical analysis, we shall perform discrete steps in rapidity according to |0J. It is thus important to 
notice that, due to the contribution at T = 0, -FTo;Ay(0) = exp[— To/(77Ay)] becomes close to 1 when Tq <C r]AY. 
As a consequence, for sufficiently small values of T, the effect of the noise will be to set the amplitude to with a 
probability almost equal to 1. This mechanism, having an effect close to a cut-off on T as applied by Brunet and 
Derrida jl^ for the sFKPP equation, ensures the compacity of the front. 

Before going on with the practical analysis, let us point out that, in order to reach very high values of the rapidity, 
we have, when necessary, worked in a comoving frame. This means that me have translated the value of Lmin and 
Lma.^ together with the propagating front. The value of Lmin has been adjusted in order to have T{L^in) > 6. The 
compactness of the front allows precise determination of the front in the ultraviolet region and the fact that the noise 
does not produce sizeable effects in the saturated region allows to handle the infrared part of the front. 



B. Simulation parameters 

We also have to discuss the discretisation of space and time. As mentioned earlier, the size of our steps in rapidity 
is related to the minimum value of the amplitude we want to describe. It is therefore important to set this to a small 
value. We have adopted a step size Ay = 0.001. The space discretisation is more complicated due to the presence 
of the Dirac distribution in L = log(A;) in the noise dispersion. For a given discretisation of size Ai, the typical 
amplitude of the noise in one step of the evolution will be proportional to ^ AY/ AL. For the optimal stability in 
numerical simulations, it is thus safe to take a value of AL small enough but still larger than the time discretisation. 
We have obtained that AL = 0.5 was a good choice. 

All our simulations have been performed with a = 0.2 which is typical for real-life situations and we have varied 
the strength aj.K of the fiuctuation term by simulating different values of k. The studies at fixed coupling are relevant 
in this approach since there is no unique way to introduce running coupling effects at the leading order. In addition, 
by fixing as, one could directly compare our results with the analytical predictions. For the purposes of numerical 
studies, we have worked with k — ^§tK instead of k. 

The last thing that needs to be fixed is the initial condition. For all the following analysis, we have chosen a 
physically-motivated amplitude T(Y = 0) which already exhibits saturation in the infrared region: 

TiY^O,k)^H-^^ (14) 
I 1 — -J otherwise. 



C. Results 



We have concentrated our studies on three major points. First, we have checked that, on an event-by-event basis, 
the dispersion of the front increases like Vy. Then, we have studied at very large rapidities the asymptotic velocity of 
the front. Finally, we have looked at the front averaged over a large number of events and we shall discuss its various 
properties. 

We be gin ou r analysis with the check that the position of the front is diffused like Vy. This diffusion is presented 
in figure 2(a) where we draw 10 events generated wit k — 1 together with their average for different values of the 
rapidity between and 50. It is obvious that, as rapidity increases, these events become more and more spread. If 
we extract the position of each front by solving T{L, Y) = 1/2 w.r.t. L for different values of Y, one can compute the 
statistical dispersion a of our set of events as a function of the rapidity. The result for is drawn in figure |2(b) for 
a sets of 10000 events with different values of k, together with a linear fit for Y > 25. This clearly proves that the 
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(a) The dashed curves show 10 events generated with 
K = 1 at y = 0, 12.5, 25, 37.5 and 50 and the plain curve 
is their average. 



(b) The points represent the square of the dispersion of the 
events as a function of rapidity together with a hnear fit at 
high rapidities. 



FIG. 2: Dispersion of the events due to the noise eflects as rapidity increases. 
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FIG. 3: Asymptotic speed of the traveling front. We observe a slow down w.r.t. the equation without fluctuations, more and 
more important as the strength of the noise is important. 



dispersion of the front position goes like &i large Y . Note that we have not gone to sufficiently low value of k to 
observe the asymptotic behaviour of the dispersion expected from Moreover, in the early stages of the evolution 
(0 < y < 10), we see that there is very few diffusion. In some sense, the diffusive process only starts to become 
important once the wavefront has reached its critical form. This point is very important since the dispersion of the 
events leads to violations of the geometric scaling. If, in the early stages of the evolution, we have very few diffusion, 
then geometric scaling will still be valid. This can thus explain why geometric scaling is indeed observed at HERA 
[T^ where Y Note however that this may of course depend on our choice of initial condition. 

As a second test, we shall check that as k increases, the asymptotic velocity of the front decreases. Since the front 
velocity has rather large pre-asymptotic corrections, it is not sufficient, as for the determination of the diffusion, to 
consider a large number of events at y « 50. The technique here is to use the comoviiig frame method to go to very 
large values of the rapidity. The speed can then be defined by 



= lim 



{L{Y)) - {L{Y,)) 

Y-Yn 



where L{Y) is the position of the front at rapidity Y, determined by T{Ls,Y) — 1/2, and Yq is a reference rapidity 
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sufRciently high for the front to be stable. PracticaUy, we have set Yq = 1000 and evolved the initial front up to 

Y = 25000. Note that here, just a few events are needed in order to determine the velocity. The results we obtained 
are presented in figure |31 for different values of the noise strength ft. This figures shows clearly the fluctuations in 
speed and the slow down w.r.t. the BK equation where the fluctuations are not taken into account. This decrease of 
the speed becomes more and more important as the constant k grows. The logarithmic dependence on k cannot really 
be obtained from these analysis since we have not been to sufficiently small values of k. This was however out of the 
scope of this paper since we simply want to check the decrease of the asymptotic speed when the noise becomes more 
and more important. Note also that the set of 10000 events used in figure 2(b) for k = 1 gives a velocity of 0.912 at 

Y = 50. Although this speed is smaller than the BK critical speed Vc, the asymptotic value ((HJ is not yet reached. 
The main motivation of this paper was to study the shape of the traveling front for physical strength of the noise 

term and its deviations from the results obtained in the absence of fluctuations. To do such an analysis, one need 
to generate a large number of events and take their average. Let us first study the evolution of the front with 
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FIG. 6: Comparison of the numerical predictions for (T) at F = 20 with the ampUtude obtained from BK evolution and with 
the sFKPP model We have used k — 1. 
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FIG. 7: Numerical results for the correlation {T^{k)) at Y ^ 40. We show (T) and (T)^ for comparison. 

rapidity. We show in figure 01 the amplitude^ (T) obtained from 10000 events up to y = 50 and different values 
of the noise parameter k. It is obvious on that figure that, as rapidity increases, the tail which is close to the BK 
result at small rapidities, starts to diverge from it due to diffusion in the events. The somehow surprising point is 
that, in the beginning of the evolution, the front including fluctuations travels faster than the BK front. This is not 
in contradiction with the fact that the asymptotic velocity needs to be smaller than the BK one, as we have seen 
previously that we need to go to much higher values of the rapidity in order to reach the asymptotic speed. The effect 
is really present in the early stages of the evolution and increases with k as we can see on the front for y = 10 and 
Y = 20. Already at y = 50 one can observe a slow down of the front which is also more important for larger values of 
K. It is also interesting to point out that analytical studies, only controlling the asymptotic behaviour (large Y and 



^ Note that the recoil observed in the last curve at the ultraviolet border of the domain is probably a non-physical numerical edge effect. 
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weak noise) , does not say anything about these early stages of the evolution. 

To conclude our analysis, let us compare our simulations with the approximation introduced by lancu and Tri- 
antafyllopoulos (see eqs. 1)10(1 and ((ll)) ). In figure El we compare the numerical results for (T) with the front 
obtained from BK evolution and with equation HlOfl where we have adjusted the parameters Qs and a in the range 
10^ < k < 10^*. We see a very good agreement between our results and the prediction (|10|l for the tail of the front. 
The description in the saturated part of the front is not the same for the simple reason that the approximation 1)10(1 
assumes that the amplitude saturates to 1 instead of taking into account the logarithmic evolution which does not 
have any consequence on the high-fc behaviour. 

Finally, we have plotted in figure [7| the average value of as a function of the momentum k ioi Y = 40 and k = 1, 
together with (T) and (T)^. In the small-fc region where the system is saturated, we clearly have {T'^) — (T)'^ which 
is coherent with the fact that we are in the dense regime and that the description is well reproduced also by the 
BK equation. In the dilute regime, we see that (T"^) becomes closer to (T) as predicted by (|ll|l . The fact that the 
adequacy is not perfect simply comes from the reason that the asymptotic regime, strictly valid for y — > oo is not 
completely reached. In that region, T being smaller than 1 with a large probability, we have < T and we expect 
(r^) to reach (T) from below. 

V. CONCLUSIONS AND DISCUSSION 

In this paper, we have studied numerically the solutions of the recently-proposed high-energy-QCD evolution equa- 
tion including fluctuations. If we perform a coarse-graining in impact-parameter and go to momentum space, this 
hierarchy becomes equivalent to a Langevin equation. To analyse the effects due to the fluctuations, we have com- 
pared the results of the complete equation with the predictions from the BK evolution equation which only takes into 
account the BFKL and saturation contributions. Instead of solving the complicated infinite hierarchy of equations, 
we have simulated the equivalent Langevin problem which allows to study the QCD amplitudes on an event-by-event 
basis ant to obtain the average amplitude by generating a large number of events. Note that we have studied the com- 
plete QCD equation, with physically acceptable values of the parameters and without performing the approximations 
leading to the analogy with statistical physics. 

Our analysis starts by showing that different events lead to a diffusion in position. More precisely, as rapidity 
increases, the positions of different events get spread with a dispersion proportional to \/(Y) and increasing with 
the importance of the noise term. In addition, it seems that this diffusion effect only becomes important once the 
wavefront is formed and hence we do not observe large dispersion in the early stages of the evolution. This may explain 
the fact that geometric scaling is indeed observed in experimental data as they lay in a region where dispersion appears 
to be small. This feature however may depend on our choice of initial condition H14|l. Indeed, to observe geometric 
scaling one has to be in a region of rapidity sufficiently large for the front to be formed and not too large so that 
fluctuations effects are not yet important. At asymptotically weak noises, the time of formation of the front is of order 
/3~^7~^ log(Ka^) with /3 « 48.2. This is faster than the equivalent time for the BK equation and one might expect 
that geometric scaling sets in early. However, this is only known for asymptotically weak noises while the effect of 
fluctuations is only known at asymptotically high rapidities. Therefore, very few things can be predicted concerning 
these competitive effects and complementary studies are required. 

We have then studied the asymptotic speed of the traveling front i.e. the rapidity evolution of the saturation scale. 
This requires to go to very large values of the rapidity in order not to be sensitive to the sub-asymptotic corrections 
and to the diffusion due to the noise which introduce corrections proportional to Y~^^^. As predicted by the sFKPP 
equation, the asymptotic saturation exponent is significantly smaller than the one obtained from the BK equation. 
Note however that since we have not studied asymptotically small values of the noise strength parameter a^K which is 
out of the scope of this paper, we have not observed the log~^(a^/t) behaviour of the corrections. It appears however 
that the slow down of the front predicted from statistical physics still holds in the QCD case. 

Finally, we have studied the interesting problem of the shape of the travelling front in the presence of fluctuations. 
This is performed by considering a large number of events and comparing the results, with different values of the 
noise parameter k, with the front obtained from BK evolution. We have observed at large rapidities the violations of 
geometric scaling in the leading part of the tail and, at very large momentum fc, a strong damping in agreement with 
the compacity property of the front in such type of stochastic evolution equation. We have observed that, in the early 
stage of the evolution the front goes faster than the front obtained from BK evolution. At higher rapidities, we already 
observe that the speed of the average front goes slower than the front without noise term. This speed also turns out 
to be smaller for larger contributions from fluctuations. Moreover, we have checked that the front predicted by the 
stochastic QCD evolution shows a behaviour consistent with the form predicted by lancu and Triantafyllopoulos. At 
high rapidities, we have checked that, according to the theoretical predictions, the average values (T^) behaves like 

(T)^ in the saturated region and converges to (T) in the tail of the traveling front. 
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In the future, it would be interesting to see if these numerical studies can help to distinguish some effects due to 
fluctuations in the experimental data. To perform these phenomenological studies, one has to consider the effects 
of running coupling. One may notice that the compacity property of the traveling front gives an efficient way to 
introduce running coupling in the evolution equation. Indeed, the dynamics of the evolution is mainly located in 
the region of the front and thus one should use as the scale for the running coupling. As for the case of the BK 
equation, one expects a front evolving like \/Y instead of Y. The diffusion of the events as a function of rapidity is 
unclear and requires more precise analyses. 

Then, one may also introduce the impact-parameter dependence which is known to lead to a stochastic equation 
with a multi-variable noise, poorly known in theoretical grounds. However, if one consider the BK equation, the 
introduction of the impact-parameter dependence [l£!| does not modify the asymptotic properties observed in the 
6-independent case. Thus we do not expect that the full hierarchy will lead to huge qualitative modifications of the 
present picture. 

At the level of purely theoretical studies, our results suggest to analyse the link with the DLA regime which appears 
to be closely related with the compacity properties of the front. 
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APPENDIX A: EQUIVALENCE BETWEEN THE LANGEVIN EQUATION AND THE HIERARCHY 

Instead of considering the full QCD case, we shall show how a Langevin equation translates into an infinite hierarchy 
of equation of the same form as in a simplified model. Since the introduction of space variable with sophisticated 
linear kernel simply introduces complications with the notations, let us consider the 0-dimensional problem. 

Our starting point is the Langevin equation 



dtu{t) = u{t) - u^{t) + y/2au{t)r](t), 

with the noise ri{t) satisfying {ri{t)) — and {ri{t)ri{t')) — S{t — t'). Since a Langevin process is not differentiable, one 
has to understand this equation as the continuum limit with an Ito prescription. If we discretise the time in intervals 
of size e, it gives 

= Uj — Uj + ^J2aUj■qj, 



where the noise satisfying 



(77j)=0 and {rijVi') = \Sjj'■ 



e 



Let us now consider the time evolution of the average value {F{u)) of a function F(u). Using again the Ito prescription 
and expanding to second order around Uj, it is straightforward to get 

dt{F{u)) = (uF'iu)) - (u^F'iu)) +a(uF"{u)). 

Taking F{u) = m" and defining m^*^-* — (u'^) we obtain the infinite hierarchy of equations 

ftu^*^) = fcu^*^) - fcuC^+i) + ak{k ~ l)u'^^-^\ 

which reproduces a linear, a merging and a splitting term. 
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